Priority Effects

# PROLOG   ##############################################################################

# PROJECT: Priority Effects Experiment 
# PURPOSE: Determine how species arrival via seed additions affects the community assembly of restored prairie plant communities
#
# AUTHOR:  Katherine Wynne (Wynnekat@msu.edu)
# COLLAB:  L. Sullivan
# CREATED: 06 December 2022
#
# FILES:   1) inner_cover_summer_2022.xlsx - vegetative cover taken in the Summer (end of June)
#          2) inner_cover_fall _2022.xlsx - vegetative cover taken at peak biomass (end ofAugust)
#          3) Checklist of Missouri Flora.xlsx - contains the information about all plants found in MO
#          4) Cleaned_PE_Detailed_Species_List_2022.xlsx - dataset containing the information about only                                                                plants found in the study
#         
#LAST UPDATED: 06 September 2023
# PROLOG   ##############################################################################

Treatment information & Setup

NON = None (no species seeded)

SIM = Simulataneous (36 species seeded on March 22nd, 2021)

LE = Lumped early (18 early-dispersing species seeded on March 23rd, 2021 then 18 late-dispersing species seeded on September 5th, 2021)

LL = Lumped late (18 late-dispersing species seeded on March 23rd, 2021 then 18 early-dispersing species seeded on September 5th, 2021)

NAT = Natural; species added according to peak dispersal date

    - May 24th 2021 (VIOPED, PACPLA)
    - June 4th 2021 (SISCAM, SPHOBT, CORLAN)
    - June 22nd 2021 (LOBSPI, TRAOHI, CARBUS, HEURIC)
    - July 11th 2021 (CORPAL, DODMEA)
    - July 25th 2021 (KOEMAC, AMOCAN)
    - August 2nd 2021 (LINSUL, MELVIR)
    - August 21st 2021 (DALCAN, DALPUR, ACHMIL)
    - September 17th 2021 (CROSAG)
    - October 4th 2021 (CHAFAS, MONFIS, RATPIN, SPOHET, RUDHIR)
    - October 18th 2021 (BIDARI, HELMOL, LESCAP)
    - November 1st 2021 (PENDIG, ERYYUC, HYPPUN, SORNUT, SCHSCO, LIAPYC)
    - November 19th 2021 (CORTRI, PYCTEN, SOLRIG)
  

Load libraries

Import Datasets

Dataset cleaning

2021

### Remove unnecessary columns

#### Get rid of the notes, unknown, and senensced columns

### Fall
inner_cover_fall_2021 <- inner_cover_fall_2021[,-c(6,7,8)]


#### change log
 # - senesced PLASPP changed to PLAVIR since PLALAN and PLAMAJ do not senesce until Oct.  (Sep 6, 2023)
 # - JUNSPP changed to JUNINT since JUNINT has been observed in 2022 and 2023 (Sep 6, 2023)


### Get rid of EMP

### Fall
inner_cover_fall_2021 <- subset(inner_cover_fall_2021, Treatment != "EMP")

2022

2023

# Make a unique identifier for year and season

### 2021

#### Season
inner_cover_fall_2021$Season <- rep("Fall", nrow=(inner_cover_fall_2021))
#### Year
inner_cover_fall_2021$Year <- rep("2021", nrow=(inner_cover_fall_2021))


### 2022

# Already has the identifiers
inner_cover_fall_2022$Year <- as.factor(inner_cover_fall_2022$Year)
inner_cover_summer_2022$Year <- as.factor(inner_cover_summer_2022$Year)

### 2023

#### Season
inner_cover_summer_2023$Season <- rep("Summer", nrow=(inner_cover_summer_2023))
inner_cover_fall_2023$Season <- rep("Fall", nrow=(inner_cover_fall_2023))

#### Year
inner_cover_fall_2023$Year <- rep("2023", nrow=(inner_cover_fall_2023))
inner_cover_fall_2023$Year <- as.factor(inner_cover_fall_2023$Year)

inner_cover_summer_2023$Year <- rep("2023", nrow=(inner_cover_summer_2023))
inner_cover_summer_2023$Year <- as.factor(inner_cover_summer_2023$Year)

# Join years together

full_2021_data <- inner_cover_fall_2021

full_2022_data <- full_join(inner_cover_summer_2022,inner_cover_fall_2022)

## Joining with `by = join_by(Season, Year, Block, Treatment, SPP6, Percent_Cover,
## Scientific_Name)`

full_2023_data <- full_join(inner_cover_summer_2023,inner_cover_fall_2023)

## Joining with `by = join_by(Block, Treatment, SPP6, Percent_Cover,
## Scientific_Name, Season, Year)`

# Join all years together

### First 2021 and 2022
full_21_22_data <- full_join(full_2021_data, full_2022_data)

## Joining with `by = join_by(Block, Treatment, SPP6, Percent_Cover,
## Scientific_Name, Season, Year)`

### Then 2021 and 2022 with 2023
full_year_data <- full_join(full_21_22_data, full_2023_data)

## Joining with `by = join_by(Block, Treatment, SPP6, Percent_Cover,
## Scientific_Name, Season, Year)`

# lump certain taxa together 

### Lump the Carex for now (*****Ask Lauren about what to do about when the CARSPP > CARFRA or CARBRE?***)

## MELOFF + MELALB -> MELSPP
## LEPVIR -> LEPSPP
## CARFRA -> CARSPP
## CARBRE -> CARSPP

for(i in 1:nrow(full_year_data)) {
  if(full_year_data[i,3] == "CARFRA"){full_year_data [i,3] <- "CARSPP"}
  if(full_year_data[i,3] == "CARBRE"){full_year_data [i,3] <- "CARSPP"}
  if(full_year_data[i,3] == "MELOFF"){full_year_data [i,3] <- "MELSPP"}
  if(full_year_data[i,3] == "MELALB"){full_year_data [i,3] <- "MELSPP"}
  if(full_year_data[i,3] == "LEPVIR"){full_year_data [i,3] <- "LEPSPP"}
  
  # Make all cover data that was less than 1 = to 1 (to make the data discrete)
  if(full_year_data[i,4] < 1) {full_year_data [i,4] <- 1}
  
}

inner_cover_max <- full_year_data  %>%
      group_by(Block, Treatment, Year, SPP6) %>%
      summarise(max_cover=max(Percent_Cover))

## `summarise()` has grouped output by 'Block', 'Treatment', 'Year'. You can
## override using the `.groups` argument.

##### Note I manually checked species that were found in both the datasets (ACHMIL, TRIREP, CORLAN) and confirmed that the code above took the highest cover value for a species seen twice


### Remove Bare
inner_cover_max_only <- subset(inner_cover_max, SPP6 != "BARE")
### Remove Litter
inner_cover_max_only  <- subset(inner_cover_max_only, SPP6 != "LITTER")

head(inner_cover_max_only) 

## # A tibble: 6 × 5
## # Groups:   Block, Treatment, Year [1]
##   Block Treatment Year  SPP6   max_cover
##   <dbl> <chr>     <chr> <chr>      <dbl>
## 1     1 LE        2021  ACAVIR         1
## 2     1 LE        2021  ACHMIL         1
## 3     1 LE        2021  AMATUB         1
## 4     1 LE        2021  AMBART        20
## 5     1 LE        2021  BARVUL         1
## 6     1 LE        2021  BROJAP         1

### Bare dataset

inner_bare_cover <- subset(inner_cover_max, SPP6 == "BARE")

### Litter dataset

inner_litter_cover <- subset(inner_cover_max, SPP6 == "LITTER")


### Manipulate 2021 dataset to match 2022 and 2023

#### Bare Ground

bare_2021 <- subset(bare_cover_2021, Cover_Type != "Litter")
bare_2021 <- subset(bare_2021, Treatment != "EMP")
names(bare_2021) <- c("Block", "Treatment", "SPP6", "max_cover")
bare_2021$Year <- rep("2021", nrow(bare_2021))

for(i in 1:nrow(bare_2021)) {
  if(bare_2021[i,3] == "Bare"){bare_2021[i,3] <- "BARE"}
}

#### Litter 

litter_2021 <- subset(bare_cover_2021 , Cover_Type != "Bare")
litter_2021 <- subset(litter_2021, Treatment != "EMP")
names(litter_2021) <- c("Block", "Treatment", "SPP6", "max_cover")
litter_2021$Year <- rep("2021", nrow(litter_2021))


for(i in 1:nrow(litter_2021)) {
  if(bare_2021[i,3] == "Litter"){bare_2021[i,3] <- "LITTER"}
}


### Full_join datasets

inner_bare_cover <- full_join(inner_bare_cover, bare_2021)

## Joining with `by = join_by(Block, Treatment, Year, SPP6, max_cover)`

inner_litter_cover <- full_join(inner_litter_cover, litter_2021)

## Joining with `by = join_by(Block, Treatment, Year, SPP6, max_cover)`

### Add species information to the dataset

#### N/A is N = native, A = non-native, G = Genus

inner_cover_max_only <- full_join(inner_cover_max_only, species_list_PE)

## Joining with `by = join_by(SPP6)`

#Remove unnecessary columns
inner_cover_max_reduced <- inner_cover_max_only[, -c(6,7,10)]

Exploratory Analysis of Diversity

Summary Table of Diversity Indices

Below is a summary table showing the mean and SD for mean C value (C), species richness (SR), and Shannon diversity index (SDI) for each seeding treatment.

Analyzing diversity

Mean C

Mean C - Model assumptions

Distribution of C values for all treatments were right skewed, indicating the majority of native plants in all treatments were of low conservatism value. SIM had the most “even” distribution of C values and almost every value was represented in this treatment.

## Warning in geom_density(bins = 30, alpha = 0.4): Ignoring unknown parameters:
## `bins`

# A couple of outliers in the NAT treatment but likely alright since the variablity in outcomes seems to be a characteristic of NAT

# only one extreme value in LE 2021 (likely because of DALCAN and VIOPED)
Summary_table %>%
  group_by(Treatment, Year) %>%
  identify_outliers(Mean_C)

## # A tibble: 4 × 8
##   Treatment Year  Block Species_richness Shannon Mean_C is.outlier is.extreme
##   <chr>     <chr> <dbl>            <int>   <dbl>  <dbl> <lgl>      <lgl>     
## 1 NAT       2023      4               18    2.44   2.69 TRUE       FALSE     
## 2 NON       2023      6               19    1.99   1.18 TRUE       FALSE     
## 3 SIM       2022      5               20    1.83   3.13 TRUE       FALSE     
## 4 SIM       2023      2               27    2.41   3.42 TRUE       FALSE

# Mean C fits a normal distribution

Summary_table  %>%
  group_by(Treatment, Year) %>%
  shapiro_test(Mean_C)

## # A tibble: 10 × 5
##    Treatment Year  variable statistic      p
##    <chr>     <chr> <chr>        <dbl>  <dbl>
##  1 LE        2022  Mean_C       0.926 0.548 
##  2 LE        2023  Mean_C       0.972 0.907 
##  3 LL        2022  Mean_C       0.890 0.317 
##  4 LL        2023  Mean_C       0.861 0.193 
##  5 NAT       2022  Mean_C       0.901 0.379 
##  6 NAT       2023  Mean_C       0.935 0.616 
##  7 NON       2022  Mean_C       0.928 0.566 
##  8 NON       2023  Mean_C       0.882 0.280 
##  9 SIM       2022  Mean_C       0.816 0.0809
## 10 SIM       2023  Mean_C       0.879 0.263

ggplot(data = Summary_table , aes(x= Mean_C)) +geom_histogram(bins =15) 

# Vast majority of points fall along the reference line

ggqqplot(Summary_table , "Mean_C", ggtheme = theme_bw()) +
  facet_grid(Year~ Treatment, labeller = "label_both")

Mean C - Model fitting

Mean C values were significantly higher in all treatments compared to the NON treatment.

#Dropped the mixed model due to singular fit -> the variance of the random effect block is estimated close to zero and does not further inform the data. 

Summary_table$Block <- as.factor(Summary_table$Block)

# Normal interaction model

#Mean_C.mod.lmer.interaction <- lmer(Mean_C~Treatment*Year+(1|Block), data =Summary_table)
#AIC(Mean_C.mod.lmer.interaction)


# Normal - additive model (has the lowest AIC)


Mean_C.mod.lm.additive <- lm(Mean_C~Treatment+Year, data =Summary_table)

## Summary results
summary(Mean_C.mod.lm.additive)

## 
## Call:
## lm(formula = Mean_C ~ Treatment + Year, data = Summary_table)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.99511 -0.34068 -0.07964  0.27121  1.09720 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    1.6743     0.1583  10.575 9.11e-15 ***
## TreatmentLL    0.4931     0.2044   2.412 0.019270 *  
## TreatmentNAT  -0.3604     0.2044  -1.763 0.083548 .  
## TreatmentNON  -1.3663     0.2044  -6.684 1.34e-08 ***
## TreatmentSIM   0.7765     0.2044   3.799 0.000371 ***
## Year2023       0.2812     0.1293   2.175 0.034032 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5007 on 54 degrees of freedom
## Multiple R-squared:  0.7199, Adjusted R-squared:  0.694 
## F-statistic: 27.76 on 5 and 54 DF,  p-value: 8.389e-14

## AIC test
AIC(Mean_C.mod.lm.additive)

## [1] 94.93578

## ANOVA type 3
anova(Mean_C.mod.lm.additive)

## Analysis of Variance Table
## 
## Response: Mean_C
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## Treatment  4 33.612  8.4029 33.5209 4.716e-14 ***
## Year       1  1.186  1.1858  4.7304   0.03403 *  
## Residuals 54 13.537  0.2507                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## I removed the random effect block because it had no effect on the model fit and explained no additional variance. (chose the most parsimonius model)

# r.squaredGLMM(Mean_C.mod.lmer.additive)


## Inclusion of the random effect marginally increased model fit 



est.lm.treatment <- emmeans::emmeans(Mean_C.mod.lm.additive, "Treatment")


#Posthoc test

## Treatment
pairs(est.lm.treatment, adjust = "tukey")

##  contrast  estimate    SE df t.ratio p.value
##  LE - LL     -0.493 0.204 54  -2.412  0.1276
##  LE - NAT     0.360 0.204 54   1.763  0.4054
##  LE - NON     1.366 0.204 54   6.684  <.0001
##  LE - SIM    -0.777 0.204 54  -3.799  0.0033
##  LL - NAT     0.853 0.204 54   4.176  0.0010
##  LL - NON     1.859 0.204 54   9.097  <.0001
##  LL - SIM    -0.283 0.204 54  -1.387  0.6388
##  NAT - NON    1.006 0.204 54   4.921  0.0001
##  NAT - SIM   -1.137 0.204 54  -5.562  <.0001
##  NON - SIM   -2.143 0.204 54 -10.483  <.0001
## 
## Results are averaged over the levels of: Year 
## P value adjustment: tukey method for comparing a family of 5 estimates

# Residual plot looks decent
plot(Mean_C.mod.lm.additive)

Mean C - Figures

Summary_table_mean$Treatment <- factor(Summary_table_mean$Treatment, levels=c("NON", "NAT", "LE", "LL", "SIM"))

lower.c <- Summary_table_mean$C_Mean - Summary_table_mean$C_SD
upper.c <- Summary_table_mean$C_Mean + Summary_table_mean$C_SD

pd <- position_dodge(width = 0.8)

C.plot <- ggplot(data = Summary_table_mean, aes(x = Year, y = C_Mean, color = Treatment))+
  geom_point(data = Summary_table_mean, size = 3, aes(shape = Treatment, fill = Treatment), position = pd)+
  geom_errorbar(aes(ymin = lower.c , ymax = upper.c), width = 0.5, position = pd) +
  labs(y = "Mean C", x = "Year" )+
  theme_classic() +
  theme(text=element_text(size=16), legend.key.size=unit(1, "cm"), legend.position="none")+
  scale_color_manual(name = "Treatment", values=c("#E69F00", "#F0E442",  "#009E73", "#56B4E9", "#0072B2"),   breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
  scale_fill_manual(name = "Treatment", values=c("#E69F00", "#F0E442",  "#009E73", "#56B4E9", "#0072B2"),   breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
   scale_shape_manual(name = "Treatment", values=c(21, 22,23, 24, 25),   breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))

C.plot

Species Richness

SIM treatments had the greatest species richness compared to all other treatments, while NON and NAT had the lowest.

SR - Model assumptions

# A couple of outliers in the NAT treatment but likely alright since the variablity in outcomes seems to be a characteristic of NAT

Summary_table %>%
  group_by(Treatment, Year) %>%
  identify_outliers(Species_richness)

## # A tibble: 5 × 8
##   Treatment Year  Block Species_richness Shannon Mean_C is.outlier is.extreme
##   <chr>     <chr> <fct>            <int>   <dbl>  <dbl> <lgl>      <lgl>     
## 1 LE        2023  3                   15    2.10   2.2  TRUE       FALSE     
## 2 LL        2022  2                   25    2.31   2.06 TRUE       FALSE     
## 3 NAT       2022  4                   22    2.18   2.31 TRUE       TRUE      
## 4 NAT       2022  6                   15    2.10   1.09 TRUE       FALSE     
## 5 NON       2023  6                   19    1.99   1.18 TRUE       FALSE

# Species richness fits a normal distribution

Summary_table %>%
  group_by(Treatment, Year) %>%
  shapiro_test(Species_richness)

## # A tibble: 10 × 5
##    Treatment Year  variable         statistic      p
##    <chr>     <chr> <chr>                <dbl>  <dbl>
##  1 LE        2022  Species_richness     0.996 0.998 
##  2 LE        2023  Species_richness     0.847 0.149 
##  3 LL        2022  Species_richness     0.814 0.0778
##  4 LL        2023  Species_richness     0.825 0.0969
##  5 NAT       2022  Species_richness     0.836 0.122 
##  6 NAT       2023  Species_richness     0.884 0.286 
##  7 NON       2022  Species_richness     0.960 0.817 
##  8 NON       2023  Species_richness     0.955 0.784 
##  9 SIM       2022  Species_richness     0.893 0.332 
## 10 SIM       2023  Species_richness     0.873 0.238

ggplot(data = Summary_table, aes(x= Species_richness)) +geom_histogram(bins =15) 

# Vast majority of points fall along the reference line

ggqqplot(Summary_table, "Species_richness", ggtheme = theme_bw()) +
  facet_grid(Year~ Treatment, labeller = "label_both")

SR - Model fitting

# Model Fitting

## Despite being count data, the distribution of Species richness is remarkabley balanced and not skewed. Model comparison using AIC reveals that the Normal model works better than the Poisson model. Model results also differ likely because of under dispersion. The more complicated Quasipoisson model that accounts for under dispersion produces similar results to the Normal mixed model. Therefore, I'm going to go with the less complicated Normal model.


SR.mod.normal <- lmer(Species_richness~Treatment+Year+(1|Block), data =Summary_table)
summary(SR.mod.normal)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Species_richness ~ Treatment + Year + (1 | Block)
##    Data: Summary_table
## 
## REML criterion at convergence: 283.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3018 -0.7002  0.0875  0.6078  1.6946 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Block    (Intercept) 2.378    1.542   
##  Residual             7.355    2.712   
## Number of obs: 60, groups:  Block, 6
## 
## Fixed effects:
##              Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   19.9333     1.0639 20.8212  18.736 1.64e-14 ***
## TreatmentLL    0.3333     1.1072 49.0000   0.301 0.764634    
## TreatmentNAT  -1.7500     1.1072 49.0000  -1.581 0.120398    
## TreatmentNON  -4.3333     1.1072 49.0000  -3.914 0.000280 ***
## TreatmentSIM   4.2500     1.1072 49.0000   3.839 0.000355 ***
## Year2023      -1.7000     0.7002 49.0000  -2.428 0.018911 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmLL TrtNAT TrtNON TrtSIM
## TreatmentLL -0.520                            
## TreatmntNAT -0.520  0.500                     
## TreatmntNON -0.520  0.500  0.500              
## TreatmntSIM -0.520  0.500  0.500  0.500       
## Year2023    -0.329  0.000  0.000  0.000  0.000

## Treatment and Year significant
## Interaction between Treatment and Year not significant

anova(SR.mod.normal, type = 3)

## Type III Analysis of Variance Table with Satterthwaite's method
##           Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## Treatment 474.77  118.69     4    49 16.1381 1.704e-08 ***
## Year       43.35   43.35     1    49  5.8941   0.01891 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Marginal R^2 (variance explained by ONLY fixed effects) = 0.4547589
## Conditional R^2 (variance explained by the entire model) = 0.5554578

## Inclusion of the random effect increased model fit 


r.squaredGLMM(SR.mod.normal)

## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.

##            R2m       R2c
## [1,] 0.4743091 0.6027586

#Posthoc test

est.lmer.SR.treatment <- emmeans::emmeans(SR.mod.normal, "Treatment")


## Treatment
pairs(est.lmer.SR.treatment , adjust = "tukey")

##  contrast  estimate   SE df t.ratio p.value
##  LE - LL     -0.333 1.11 49  -0.301  0.9981
##  LE - NAT     1.750 1.11 49   1.581  0.5166
##  LE - NON     4.333 1.11 49   3.914  0.0025
##  LE - SIM    -4.250 1.11 49  -3.839  0.0031
##  LL - NAT     2.083 1.11 49   1.882  0.3406
##  LL - NON     4.667 1.11 49   4.215  0.0010
##  LL - SIM    -3.917 1.11 49  -3.538  0.0076
##  NAT - NON    2.583 1.11 49   2.333  0.1519
##  NAT - SIM   -6.000 1.11 49  -5.419  <.0001
##  NON - SIM   -8.583 1.11 49  -7.753  <.0001
## 
## Results are averaged over the levels of: Year 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 5 estimates

# Poisson Model

# SR.mod.pois <- glmer(Species_richness~Treatment+Year+(1|Block),family = "poisson", data =Summary_table)
# summary(SR.mod.pois)
# Anova(SR.mod.pois, type = 3)

# Check dispersion

## Underdispersed! 0.659 
## Interaction between treatment and year are not significant

# Dispersion_glmer(SR.mod.pois)


# Not accounting for the underdispersion mattered a lot! .

## Used a quassipoisson model instead and year is significant

#m6 <-glmmPQL(Species_richness~Treatment+Year,
             # random = ~ 1 | Block,
             # family = quasipoisson(link='log'), 
             # data = Summary_table)

# summary(m6)
# Anova(m6, type = 3)

# Residual plot looks decent except for one large outlier
plot(SR.mod.normal)

SR - Figures

Summary_table_mean$Treatment <- factor(Summary_table_mean$Treatment, levels=c("NON", "NAT", "LE", "LL", "SIM"))

lower.SR.c <- Summary_table_mean$SR_Mean - Summary_table_mean$SR_SD
upper.SR.c <- Summary_table_mean$SR_Mean + Summary_table_mean$SR_SD

pd <- position_dodge(width = 0.8)

SR.plot1 <- ggplot(data = Summary_table_mean, aes(x = Year, y = SR_Mean, color = Treatment))+
  geom_point(data = Summary_table_mean, size = 3, aes(shape = Treatment, fill = Treatment), position = pd)+
  geom_errorbar(aes(ymin = lower.SR.c , ymax = upper.SR.c), width = 0.5, position = pd) +
  labs(y = "Total Species Richness", x = "Year" )+
  theme_classic() +
  ylim(10,30)+
  theme(text=element_text(size=16), legend.key.size=unit(1, "cm"))+
  scale_color_manual(name = "Treatment", values=c("#E69F00", "#F0E442",  "#009E73", "#56B4E9", "#0072B2"),   breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
  scale_fill_manual(name = "Treatment", values=c("#E69F00", "#F0E442",  "#009E73", "#56B4E9", "#0072B2"),   breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
   scale_shape_manual(name = "Treatment", values=c(21, 22,23, 24, 25),   breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
    theme(legend.position="bottom")

SR.plot <- SR.plot1 +
    theme(legend.position="none")


leg <- get_legend(SR.plot1)

Shannon Diversity Index

SDI - Model assumptions

SIM and LL treatments had the greatest Shannon diversity index value while NON had the lowest. NAT and LE appear comparable.

# Two extreme outliers (one in NAT 2021 and another in SIM 2021)
Summary_table %>%
  group_by(Treatment, Year) %>%
  identify_outliers(Shannon)

## # A tibble: 1 × 8
##   Treatment Year  Block Species_richness Shannon Mean_C is.outlier is.extreme
##   <chr>     <chr> <fct>            <int>   <dbl>  <dbl> <lgl>      <lgl>     
## 1 NON       2022  3                   11    1.55  0.429 TRUE       FALSE

# Except for one point, Shannon diversity fits a normal distribution

Summary_table %>%
  group_by(Treatment, Year) %>%
  shapiro_test(Shannon)

## # A tibble: 10 × 5
##    Treatment Year  variable statistic     p
##    <chr>     <chr> <chr>        <dbl> <dbl>
##  1 LE        2022  Shannon      0.971 0.901
##  2 LE        2023  Shannon      0.967 0.873
##  3 LL        2022  Shannon      0.977 0.936
##  4 LL        2023  Shannon      0.854 0.169
##  5 NAT       2022  Shannon      0.950 0.743
##  6 NAT       2023  Shannon      0.931 0.586
##  7 NON       2022  Shannon      0.848 0.151
##  8 NON       2023  Shannon      0.950 0.739
##  9 SIM       2022  Shannon      0.872 0.235
## 10 SIM       2023  Shannon      0.940 0.663

ggplot(data = Summary_table, aes(x= Shannon)) +geom_histogram(bins =15) 

# Vast majority of points fall along the reference line

ggqqplot(Summary_table, "Shannon", ggtheme = theme_bw()) +
  facet_grid(Year~ Treatment, labeller = "label_both")

SDI - Model fitting

Shannon.lmer.mod <- lmer(Shannon~Treatment+Year+(1|Block), data =Summary_table)
summary(Shannon.lmer.mod)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shannon ~ Treatment + Year + (1 | Block)
##    Data: Summary_table
## 
## REML criterion at convergence: 2.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.60837 -0.47607  0.08046  0.51127  2.04855 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Block    (Intercept) 0.007995 0.08942 
##  Residual             0.042202 0.20543 
## Number of obs: 60, groups:  Block, 6
## 
## Fixed effects:
##              Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)   2.16123    0.07452 28.51199  29.003  < 2e-16 ***
## TreatmentLL   0.05159    0.08387 49.00000   0.615  0.54133    
## TreatmentNAT -0.12219    0.08387 49.00000  -1.457  0.15152    
## TreatmentNON -0.34644    0.08387 49.00000  -4.131  0.00014 ***
## TreatmentSIM  0.12917    0.08387 49.00000   1.540  0.12994    
## Year2023     -0.01536    0.05304 49.00000  -0.290  0.77333    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmLL TrtNAT TrtNON TrtSIM
## TreatmentLL -0.563                            
## TreatmntNAT -0.563  0.500                     
## TreatmntNON -0.563  0.500  0.500              
## TreatmntSIM -0.563  0.500  0.500  0.500       
## Year2023    -0.356  0.000  0.000  0.000  0.000

anova(Shannon.lmer.mod, type = 3)

## Type III Analysis of Variance Table with Satterthwaite's method
##            Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## Treatment 1.65266 0.41316     4    49  9.7902 6.691e-06 ***
## Year      0.00354 0.00354     1    49  0.0839    0.7733    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# No significant interaction between treatment and year



## Marginal R^2 (variance explained by ONLY fixed effects) = 0.330704
## Conditional R^2 (variance explained by the entire model) = 0.3844412

## Inclusion of the random effect increased model fit 

r.squaredGLMM(Shannon.lmer.mod)

##           R2m       R2c
## [1,] 0.358652 0.4608049

# Post-Hoc comparison

pairs(emmeans::emmeans(Shannon.lmer.mod, "Treatment"))

##  contrast  estimate     SE df t.ratio p.value
##  LE - LL    -0.0516 0.0839 49  -0.615  0.9720
##  LE - NAT    0.1222 0.0839 49   1.457  0.5947
##  LE - NON    0.3464 0.0839 49   4.131  0.0013
##  LE - SIM   -0.1292 0.0839 49  -1.540  0.5420
##  LL - NAT    0.1738 0.0839 49   2.072  0.2486
##  LL - NON    0.3980 0.0839 49   4.746  0.0002
##  LL - SIM   -0.0776 0.0839 49  -0.925  0.8858
##  NAT - NON   0.2242 0.0839 49   2.674  0.0727
##  NAT - SIM  -0.2514 0.0839 49  -2.997  0.0331
##  NON - SIM  -0.4756 0.0839 49  -5.671  <.0001
## 
## Results are averaged over the levels of: Year 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 5 estimates

# A little heteroscedastic (bowtie shaped but probably ok)
plot(Shannon.lmer.mod)

SDI - Figures

Summary_table_mean$Treatment <- factor(Summary_table_mean$Treatment, levels=c("NON", "NAT", "LE", "LL", "SIM"))

lower.SDI.c <- Summary_table_mean$SDI_Mean - Summary_table_mean$SDI_SD
upper.SDI.c <- Summary_table_mean$SDI_Mean + Summary_table_mean$SDI_SD

pd <- position_dodge(width = 0.8)

SDI.plot <- ggplot(data = Summary_table_mean, aes(x = Year, y = SDI_Mean, color = Treatment))+
  geom_point(data =Summary_table_mean, size = 3, aes(shape = Treatment, fill = Treatment), position = pd)+
  geom_errorbar(aes(ymin = lower.SDI.c , ymax = upper.SDI.c), width = 0.5, position = pd) +
  labs(y = "SDI", x = "Year" )+
  theme_classic() +
  theme(text=element_text(size=16), legend.key.size=unit(1, "cm"), legend.position="none")+
  scale_color_manual(name = "Treatment", values=c("#E69F00", "#F0E442",  "#009E73", "#56B4E9", "#0072B2"),   breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
  scale_fill_manual(name = "Treatment", values=c("#E69F00", "#F0E442",  "#009E73", "#56B4E9", "#0072B2"),   breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
   scale_shape_manual(name = "Treatment", values=c(21, 22,23, 24, 25),   breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))

SDI.plot 

Plot Panel

Comparing seeded cover (in the inner 1 m^2 sampling area) for each treatment

note some of the species used in the seed mix were present in the seed bank at my study site (e.g., RUDHIR and PENDIG). Additionally, some spillover occurred with seeded species seeding into other treatments like NON (e.g., BIDARI)

Total Seeded Cover

TSC - Model fitting

## 
## Call:
## glm(formula = cbind(seeded_cover, tot_cover) ~ Treatment * Year, 
##     family = binomial, data = Total_Cover_Tot)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -6.4760  -2.0224  -0.4559   1.5998   4.5960  
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           -2.11998    0.08881 -23.870  < 2e-16 ***
## TreatmentLL            0.94324    0.10412   9.059  < 2e-16 ***
## TreatmentNAT          -0.38112    0.13799  -2.762  0.00575 ** 
## TreatmentSIM           1.06270    0.10397  10.221  < 2e-16 ***
## Year2023               0.62783    0.11282   5.565 2.63e-08 ***
## TreatmentLL:Year2023  -0.41302    0.13685  -3.018  0.00254 ** 
## TreatmentNAT:Year2023 -0.36533    0.18174  -2.010  0.04441 *  
## TreatmentSIM:Year2023 -0.37777    0.13662  -2.765  0.00569 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 837.79  on 47  degrees of freedom
## Residual deviance: 293.23  on 40  degrees of freedom
## AIC: 557.58
## 
## Number of Fisher Scoring iterations: 5

## Analysis of Deviance Table (Type III tests)
## 
## Response: cbind(seeded_cover, tot_cover)
##                LR Chisq Df Pr(>Chisq)    
## Treatment       272.445  3  < 2.2e-16 ***
## Year             32.035  1  1.514e-08 ***
## Treatment:Year   10.305  3    0.01614 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##  contrast                    odds.ratio     SE  df null z.ratio p.value
##  LE Year2022 / LL Year2022        0.389 0.0405 Inf    1  -9.059  <.0001
##  LE Year2022 / NAT Year2022       1.464 0.2020 Inf    1   2.762  0.1048
##  LE Year2022 / SIM Year2022       0.346 0.0359 Inf    1 -10.221  <.0001
##  LE Year2022 / LE Year2023        0.534 0.0602 Inf    1  -5.565  <.0001
##  LE Year2022 / LL Year2023        0.314 0.0328 Inf    1 -11.076  <.0001
##  LE Year2022 / NAT Year2023       1.126 0.1469 Inf    1   0.909  0.9853
##  LE Year2022 / SIM Year2023       0.269 0.0281 Inf    1 -12.574  <.0001
##  LL Year2022 / NAT Year2022       3.760 0.4466 Inf    1  11.150  <.0001
##  LL Year2022 / SIM Year2022       0.887 0.0680 Inf    1  -1.558  0.7751
##  LL Year2022 / LE Year2023        1.371 0.1210 Inf    1   3.572  0.0085
##  LL Year2022 / LL Year2023        0.807 0.0625 Inf    1  -2.774  0.1016
##  LL Year2022 / NAT Year2023       2.892 0.3181 Inf    1   9.654  <.0001
##  LL Year2022 / SIM Year2023       0.691 0.0534 Inf    1  -4.784  <.0001
##  NAT Year2022 / SIM Year2022      0.236 0.0280 Inf    1 -12.169  <.0001
##  NAT Year2022 / LE Year2023       0.365 0.0461 Inf    1  -7.977  <.0001
##  NAT Year2022 / LL Year2023       0.215 0.0256 Inf    1 -12.917  <.0001
##  NAT Year2022 / NAT Year2023      0.769 0.1096 Inf    1  -1.842  0.5909
##  NAT Year2022 / SIM Year2023      0.184 0.0219 Inf    1 -14.231  <.0001
##  SIM Year2022 / LE Year2023       1.545 0.1361 Inf    1   4.935  <.0001
##  SIM Year2022 / LL Year2023       0.909 0.0702 Inf    1  -1.234  0.9218
##  SIM Year2022 / NAT Year2023      3.259 0.3580 Inf    1  10.753  <.0001
##  SIM Year2022 / SIM Year2023      0.779 0.0600 Inf    1  -3.246  0.0258
##  LE Year2023 / LL Year2023        0.588 0.0523 Inf    1  -5.971  <.0001
##  LE Year2023 / NAT Year2023       2.109 0.2495 Inf    1   6.312  <.0001
##  LE Year2023 / SIM Year2023       0.504 0.0447 Inf    1  -7.729  <.0001
##  LL Year2023 / NAT Year2023       3.585 0.3958 Inf    1  11.563  <.0001
##  LL Year2023 / SIM Year2023       0.857 0.0667 Inf    1  -1.988  0.4901
##  NAT Year2023 / SIM Year2023      0.239 0.0264 Inf    1 -12.981  <.0001
## 
## P value adjustment: tukey method for comparing a family of 8 estimates 
## Tests are performed on the log odds ratio scale

TSC - Figures

Early season seeded cover

ESC - Model fitting

## 
## Call:
## glm(formula = cbind(seeded_cover, tot_cover) ~ Treatment + Year, 
##     family = binomial, data = Total_Cover_Early)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.8127  -1.5620  -0.8626   0.7913   5.1753  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -2.9526     0.1056 -27.955  < 2e-16 ***
## TreatmentLL   -3.3561     0.3411  -9.838  < 2e-16 ***
## TreatmentNAT  -2.0984     0.1995 -10.519  < 2e-16 ***
## TreatmentSIM  -0.3623     0.1070  -3.387 0.000707 ***
## Year2023       1.0821     0.1106   9.780  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 734.86  on 47  degrees of freedom
## Residual deviance: 247.79  on 43  degrees of freedom
## AIC: 381.05
## 
## Number of Fisher Scoring iterations: 6

## Analysis of Deviance Table (Type III tests)
## 
## Response: cbind(seeded_cover, tot_cover)
##           LR Chisq Df Pr(>Chisq)    
## Treatment   373.50  3  < 2.2e-16 ***
## Year        106.69  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##  contrast  odds.ratio     SE  df null z.ratio p.value
##  LE / LL      28.6759 9.7820 Inf    1   9.838  <.0001
##  LE / NAT      8.1535 1.6265 Inf    1  10.519  <.0001
##  LE / SIM      1.4366 0.1537 Inf    1   3.387  0.0039
##  LL / NAT      0.2843 0.1089 Inf    1  -3.284  0.0056
##  LL / SIM      0.0501 0.0172 Inf    1  -8.706  <.0001
##  NAT / SIM     0.1762 0.0360 Inf    1  -8.502  <.0001
## 
## Results are averaged over the levels of: Year 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log odds ratio scale

ESC - Figures

Fall dispersing species only

LSC - Model fitting

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(seeded_cover, tot_cover) ~ Treatment + Year + (1 | Block)
##    Data: Total_Cover_Late
## 
##      AIC      BIC   logLik deviance df.resid 
##    605.3    616.5   -296.6    593.3       42 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.3147 -1.7714 -0.6057  1.2574  6.5425 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Block  (Intercept) 0.008129 0.09016 
## Number of obs: 48, groups:  Block, 6
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.74863    0.09364 -29.352  < 2e-16 ***
## TreatmentLL   1.59358    0.09078  17.554  < 2e-16 ***
## TreatmentNAT  0.16919    0.11161   1.516  0.12953    
## TreatmentSIM  1.55200    0.09191  16.885  < 2e-16 ***
## Year2023      0.13908    0.05056   2.751  0.00594 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmLL TrtNAT TrtSIM
## TreatmentLL -0.796                     
## TreatmntNAT -0.643  0.664              
## TreatmntSIM -0.789  0.807  0.656       
## Year2023    -0.279  0.018 -0.002  0.024

## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: cbind(seeded_cover, tot_cover)
##                Chisq Df Pr(>Chisq)    
## (Intercept) 861.5598  1    < 2e-16 ***
## Treatment   566.5093  3    < 2e-16 ***
## Year          7.5683  1    0.00594 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##  contrast  odds.ratio     SE  df null z.ratio p.value
##  LE / LL        0.203 0.0184 Inf    1 -17.554  <.0001
##  LE / NAT       0.844 0.0942 Inf    1  -1.516  0.4278
##  LE / SIM       0.212 0.0195 Inf    1 -16.885  <.0001
##  LL / NAT       4.155 0.3536 Inf    1  16.738  <.0001
##  LL / SIM       1.042 0.0592 Inf    1   0.732  0.8841
##  NAT / SIM      0.251 0.0216 Inf    1 -16.031  <.0001
## 
## Results are averaged over the levels of: Year 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log odds ratio scale

LSC - Figures

Plot Panel

Community Analysis

NMDS Ordination

Year: 2022

Year: 2023

All

PERMANOVA

PERMANOVA analysis revealed that treatment did have a significant effect on plant community composition. Treatment also explained ~ 20% of the variation in the plant community. Post-hoc analysis shows that LE, NON, and NAT treatments were not significantly different from each other but did differ in plant community composition compared to the SIM and LL seeding treatments. Additionally, the SIM and LL treatments were only marginally different (p = 0.053). Based on the NMDS above, the LL plant community appears as a subset of the SIM plant community. Despite receiving the same pool of seeded species, species added to the plots later in the growing season did not establish as well. Barriers to establishment also appeared higher for spring/summer dispersing species compared to fall dispersing species, particularly when added without priority.

## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 1000
## 
## adonis2(formula = inner.max ~ Treatment + Year, data = inner_wide, permutations = 1000, method = "bray")
##           Df SumOfSqs      R2      F   Pr(>F)    
## Treatment  4   3.1762 0.17687 3.3089 0.000999 ***
## Year       1   1.8224 0.10148 7.5940 0.000999 ***
## Residual  54  12.9588 0.72164                    
## Total     59  17.9574 1.00000                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## $parent_call
## [1] "inner.max ~ Treatment + Year , strata = Null , permutations 999"
## 
## $LE_vs_LL
##           Df SumOfSqs      R2      F Pr(>F)    
## Treatment  1   0.9839 0.14313 4.2642  0.001 ***
## Year       1   1.0450 0.15202 4.5291  0.001 ***
## Residual  21   4.8456 0.70486                  
## Total     23   6.8745 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $LE_vs_NAT
##           Df SumOfSqs      R2      F Pr(>F)    
## Treatment  1   0.4410 0.06307 1.6582  0.040 *  
## Year       1   0.9656 0.13812 3.6310  0.001 ***
## Residual  21   5.5845 0.79881                  
## Total     23   6.9911 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $LE_vs_NON
##           Df SumOfSqs      R2      F Pr(>F)    
## Treatment  1   0.6830 0.10287 2.8438  0.001 ***
## Year       1   0.9126 0.13746 3.8000  0.001 ***
## Residual  21   5.0433 0.75967                  
## Total     23   6.6389 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $LE_vs_SIM
##           Df SumOfSqs      R2      F Pr(>F)    
## Treatment  1   0.7381 0.10692 2.9548  0.001 ***
## Year       1   0.9197 0.13323 3.6819  0.001 ***
## Residual  21   5.2454 0.75986                  
## Total     23   6.9031 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $LL_vs_NAT
##           Df SumOfSqs      R2      F Pr(>F)    
## Treatment  1   0.7043 0.10632 2.8839  0.001 ***
## Year       1   0.7917 0.11951 3.2417  0.001 ***
## Residual  21   5.1287 0.77418                  
## Total     23   6.6247 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $LL_vs_NON
##           Df SumOfSqs      R2      F Pr(>F)    
## Treatment  1   1.1298 0.17501 5.2743  0.001 ***
## Year       1   0.8276 0.12820 3.8636  0.001 ***
## Residual  21   4.4985 0.69680                  
## Total     23   6.4560 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $LL_vs_SIM
##           Df SumOfSqs      R2      F Pr(>F)    
## Treatment  1   0.5752 0.09413 2.6448  0.002 ** 
## Year       1   0.9681 0.15842 4.4511  0.001 ***
## Residual  21   4.5672 0.74744                  
## Total     23   6.1105 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $NAT_vs_NON
##           Df SumOfSqs      R2      F Pr(>F)    
## Treatment  1   0.4520 0.07021 1.8176  0.033 *  
## Year       1   0.7635 0.11859 3.0701  0.001 ***
## Residual  21   5.2222 0.81119                  
## Total     23   6.4377 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $NAT_vs_SIM
##           Df SumOfSqs      R2      F Pr(>F)    
## Treatment  1   0.9454 0.13241 3.6145  0.001 ***
## Year       1   0.7019 0.09831 2.6836  0.002 ** 
## Residual  21   5.4929 0.76928                  
## Total     23   7.1403 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $NON_vs_SIM
##           Df SumOfSqs      R2      F Pr(>F)    
## Treatment  1   1.2878 0.18695 5.5308  0.001 ***
## Year       1   0.7109 0.10321 3.0533  0.001 ***
## Residual  21   4.8897 0.70984                  
## Total     23   6.8884 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## attr(,"class")
## [1] "pwadstrata" "list"

Weather data

# Precipitation



month_tot <- weather %>% 
  filter(Year != "2021") %>% 
  group_by(Month, Year) %>% 
  summarize(tot_month_PRP_mm = sum(Total_PRP_mm))

## `summarise()` has grouped output by 'Month'. You can override using the
## `.groups` argument.

year_tot <- weather %>% 
  filter(Year != "2021") %>% 
  group_by(Year) %>% 
  summarize(tot_PRP_mm = sum(Total_PRP_mm))

library(mosaic)

## Registered S3 method overwritten by 'mosaic':
##   method                           from   
##   fortify.SpatialPolygonsDataFrame ggplot2
## 
## The 'mosaic' package masks several functions from core packages in order to add 
## additional features.  The original behavior of these functions should not be affected by this.
## 
## Attaching package: 'mosaic'
## 
## The following objects are masked from 'package:rstatix':
## 
##     cor_test, prop_test, t_test
## 
## The following object is masked from 'package:lmerTest':
## 
##     rand
## 
## The following object is masked from 'package:lme4':
## 
##     factorize
## 
## The following object is masked from 'package:Matrix':
## 
##     mean
## 
## The following object is masked from 'package:cowplot':
## 
##     theme_map
## 
## The following object is masked from 'package:permute':
## 
##     shuffle
## 
## The following objects are masked from 'package:dplyr':
## 
##     count, do, tally
## 
## The following object is masked from 'package:purrr':
## 
##     cross
## 
## The following object is masked from 'package:ggplot2':
## 
##     stat
## 
## The following objects are masked from 'package:stats':
## 
##     binom.test, cor, cor.test, cov, fivenum, IQR, median, prop.test,
##     quantile, sd, t.test, var
## 
## The following objects are masked from 'package:base':
## 
##     max, mean, min, prod, range, sample, sum

favstats(year_tot$tot_PRP_mm)

##      min      Q1 median       Q3      max    mean       sd  n missing
##  665.988 848.995 934.72 1035.304 1148.842 927.735 143.0993 10       0

# Average air temperature


library(mosaic)

favstats(weather$Average_Air_Temp_C)

##    min       Q1   median       Q3      max     mean       sd    n missing
##  -21.5 4.277778 13.33333 21.55556 32.66667 12.39584 10.45545 4018       0